% Plot the power at M2 (12.421 hours) of the observatory data on to the
% predicted power map


S = dir('/home/mnair/projects/tides/*spectra_new.mat'); % Only PMTM results (variable P) is useful

for i = 1:length(S),
    robstn(i,:) = S(i).name(1:3);
end;
load /home/mnair/projects/tides/M2_results_720_1440 Hr_power;
periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];
% period_string = ['S6';'S5';'S4';'S3';'K2';'S2';'M2';'N2';'K1';'S1';'O1'];

%OK I don't need periods < 10

periods_12 = [11.967236 12  12.421 12.6583 ];
periods_12_text = [11.9 12.01  12.43 12.67 ];
period_string_12 = ['K2';'S2';'M2';'N2'];

periods_24 = [23.934472 24 25.891];
periods_24_text = [23.8 24.01 25.9];
period_string_24 = ['K1';'S1';'O1'];

% 
% diff_ls_all = 0;
% diff_ls_nigh = 0;
% diff_rob_all = 0;
% diff_rob_night = 0;

nstsn_plot = 0;

fid = fopen('/home/mnair/projects/tides/M2_results_table.txt','wt');
fprintf(fid,'STN Longitude Latitude PMTM_Noise Observed_power_average Observed_power_1cell M2_LS_ALL M2_LS_NGHT M2_RB_ALL M2 RB_NIGHT\n');

for nstn = 1:length (S);
    
      eval(['load /home/mnair/projects/tides/' S(nstn).name ]);

    
    
    median_amp = median(sqrt(P/(en-st)));
    
    if median_amp < 1000
            
            
            
            % get the m2 power from different processing methods
            m2_power_night   = spectra_data_night_ls(7)^2;
            m2_power_all = spectra_data_all_ls(7)^2;
            m2_power_all_rob = spectra_data_all_rob(7)^2;
            m2_power_night_rob = spectra_data_night_rob(7)^2;
            
            % Average the predicted power over a square
            st_lat = round((90-latitude)*4)-2;
            en_lat = round((90-latitude)*4)+2;
            st_lon = round(longitude*4)-2;
            en_lon = round(longitude*4)+2;
            
            if st_lat < 1,
                st_lat = 1;
            end;
            
            if en_lon > 1440,
                en_lon = 1440;
            end;
            
             if en_lat > 720,
                st_lat = 720;
            end;
            
            if st_lon < 1,
                st_lon = 1;
            end;
            
            pred_m2_power = mean(mean(Hr_power(st_lat:en_lat,st_lon:en_lon))); %Average over an area
            pred_m2_power_1cell = Hr_power(round((90-latitude)*4),round(longitude*4)); % Value in the cell
            
            
            % calculate  error
            nstsn_plot = nstsn_plot + 1;
            diff_ls_all(nstsn_plot) = pred_m2_power - m2_power_all ;
            diff_ls_night(nstsn_plot) = pred_m2_power - m2_power_night ;
            diff_rob_all(nstsn_plot) = pred_m2_power - m2_power_all_rob ;
            diff_rob_night(nstsn_plot) = pred_m2_power - m2_power_night_rob ;
            
            diff_ls_all_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_all ;
            diff_ls_night_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_night ;
            diff_rob_all_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_all_rob ;
            diff_rob_night_1cell(nstsn_plot) = pred_m2_power_1cell - m2_power_night_rob ;
            
            %median_amp
            
            median_amp_array(nstsn_plot) = median_amp;
            pred_power_array(nstsn_plot) = pred_m2_power ;
            m2_power_night_rob_array(nstsn_plot) = m2_power_night_rob;
            m2_power_night_ls_array(nstsn_plot) = m2_power_night;

            % relative difference
            
            %fprintf(fid,'STN Longitude Latitude Median_Power Observed_power_average Observed_power_1cell M2_LS_ALL M2_LS_NGHT M2_RB_ALL M2_RB_NIGHT\n');
            fprintf(fid,'%s %7.2f %7.2f %10.4f %10.4f %10.4f  %10.4f  %10.4f  %10.4f %10.4f \n', ...
            S(nstn).name(1:3), ...
            longitude, ...
            latitude, ...
            median_amp, ...
            pred_m2_power, ...
            pred_m2_power_1cell, ...
            m2_power_all, ...
            m2_power_night, ...
            m2_power_all_rob, ...
            m2_power_night_rob);
            

           
        end;
    end;


fclose all
